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Abstract  Construction  of  an  accurate  theory  of  orbits  about  a  precessing  and  nutating  oblate 
planet,  in  terms  of  osculating  elements  defined  in  a  frame  associated  with  the  equator  of  date, 
was  started  in  Efroimsky  and  Goldreich  (2004)  and  Efroimsky  (2004,  2005,  2006a,  b).  Here 
we  continue  this  line  of  research  by  combining  that  analytical  machinery  with  numerical 
tools.  Our  model  includes  three  factors:  the  J2  of  the  planet,  its  nonuniform  equinoctial  pre¬ 
cession  described  by  the  Colombo  formalism,  and  the  gravitational  pull  of  the  Sun.  This 
semianalytical  and  seminumerical  theory,  based  on  the  Lagrange  planetary  equations  for 
the  Keplerian  elements,  is  then  applied  to  Deimos  on  very  long  time  scales  (up  to  1  billion 
years).  In  parallel  with  the  said  semianalytical  theory  for  the  Keplerian  elements  defined 
in  the  co-precessing  equatorial  frame,  we  have  also  carried  out  a  completely  independent, 
purely  numerical,  integration  in  a  quasi-inertial  Cartesian  frame.  The  results  agree  to  within 
fractions  of  a  percent,  thus  demonstrating  the  applicability  of  our  semianalytical  model  over 
long  timescales.  Another  goal  of  this  work  was  to  make  an  independent  check  of  whether 
the  equinoctial-precession  variations  predicted  for  a  rigid  Mars  by  the  Colombo  model  could 


We  use  the  term  “precession”  in  its  general  meaning,  which  includes  any  change  of  the  instantaneous  spin 
axis.  So  generally  defined  precession  embraces  the  entire  spectrum  of  spin-axis  variations — from  the  polar 
wander  and  nutations  through  the  Chandler  wobble  through  the  equinoctial  precession. 
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have  been  sufficient  to  repel  its  moons  away  from  the  equator.  An  answer  to  this  question,  in 
combination  with  our  knowledge  of  the  current  position  of  Phobos  and  Deimos,  will  help  us 
to  understand  whether  the  Martian  obliquity  could  have  undergone  the  large  changes  ensuing 
from  the  said  model  (Ward  1973;  Touma  and  Wisdom  1993, 1994;  Laskar  and  Robutel  1993), 
or  whether  the  changes  ought  to  have  been  less  intensive  (Bills  2006;  Paige  et  al.  2007).  It  has 
turned  out  that,  for  low  initial  inclinations,  the  orbit  inclination  reckoned  from  the  precess- 
ing  equator  of  date  is  subject  only  to  small  variations.  This  is  an  extension,  to  non-uniform 
equinoctial  precession  given  by  the  Colombo  model,  of  an  old  result  obtained  by  Goldreich 
(1965)  for  the  case  of  uniform  precession  and  a  low  initial  inclination.  However,  near-polar 
initial  inclinations  may  exhibit  considerable  variations  for  up  to  d=  10  deg  in  magnitude.  This 
result  is  accentuated  when  the  obliquity  is  large.  Nevertheless,  the  analysis  confirms  that  an 
oblate  planet  can,  indeed,  afford  large  variations  of  the  equinoctial  precession  over  hundreds 
of  millions  of  years,  without  repelling  its  near- equatorial  satellites  away  from  the  equator  of 
date:  the  satellite  inclination  oscillates  but  does  not  show  a  secular  increase.  Nor  does  it  show 
secular  decrease,  a  fact  that  is  relevant  to  the  discussion  of  the  possibility  of  high-inclination 
capture  of  Phobos  and  Deimos. 

Keywords  Orbital  elements  •  Osculating  elements  •  Mars  •  Natural  satellites  • 

Natural  satellites’  orbits  •  Deimos  •  Equinoctial  precession  •  The  Goldreich  lock 

1  Introduction 

1.1  Statement  of  purpose 

The  goal  of  this  paper  is  to  explore,  by  two  very  different  methods,  inclination  variations  of  a 
solar-gravity-perturbed  satellite  orbiting  an  oblate  planet  subject  to  nonuniform  equinoctial 
precession.  This  nonuniformity  of  precession  is  caused  by  the  presence  of  the  other  planets. 
Their  gravitational  pull  entails  precession  of  the  circumsolar  orbit  of  our  planet;  this  entails 
variations  of  the  solar  torque  acting  on  it;  these  torque  variations  make  the  planet’s  equinoctial 
precession  nonuniform;  and  this  nonuniformity,  in  its  turn,  influences  the  behaviour  of  the 
planet’s  satellites.  This  influence  is  feeble,  and  we  trace  with  a  high  accuracy  whether  it  results, 
over  aeons,  in  purely  periodic  changes  in  inclination  or  can  accumulate  to  secular  changes. 

This  work  is  but  a  small  part  of  a  larger  project  whose  eventual  goal  is  to  build  up  a 
comprehensive  tool  for  computation  of  long-term  orbital  evolution  of  satellites.  Building  this 
tool,  block  by  block,  we  are  beginning  with  only  three  components — the  planet’s  oblateness, 
the  direct  pull  of  the  Sun  on  the  satellite,  and  the  planet’s  precession.  These  phenomena 
bare  a  marked  effect  on  the  evolution  of  the  orbit.  In  our  subsequent  publications,  we  shall 
incorporate  more  effects  into  our  model — the  triaxiality,  and  the  bodily  tides. 

1.2  Motivation 

One  motivation  for  this  work  stems  from  our  intention  to  carry  out  an  independent  check 
of  whether  the  equinoctial-precession  changes  predicted  for  a  rigid  Mars  by  the  Colombo 
model  could  have  been  sufficient  to  repel  its  moons  away  from  the  equator.  An  answer  to  this 
question,  in  combination  with  our  knowledge  of  the  current  position  of  Phobos  and  Deimos, 
will  help  us  to  understand  whether  the  Martian  obliquity  variations  could  indeed  have  under¬ 
gone  the  large  variations  resulting  from  the  Colombo  model,  or  whether  the  actual  variations 
ought  to  have  smaller  magnitude.  Such  a  check  is  desirable  because  the  current,  Colombo- 
model-based  theory  of  equinoctial  precession  (Ward  1973;  Touma  and  Wisdom  1993,  1994; 


Springer 


Long-term  evolution  of  orbits 


263 


Laskar  and  Robutel  1993),  incorporates  several  approximations.  First,  the  Colombo  equa¬ 
tion  is  derived  under  the  assertion  that  the  planet  is  rigid  and  is  always  in  its  principal  spin 
state,  the  angular-momentum  vector  staying  parallel  to  the  angular- velocity  one.  Second,  this 
description,  being  only  a  model,  ignores  the  possibility  of  planetary  catastrophes  that  might 
have  altered  the  planet’s  spin  mode.  Third,  this  description  ignores  that  sometimes  even  weak 
dissipation  (caused,  for  example,  by  tides)  may  be  sufficient  to  quell  chaos  and  regularise 
the  motion,  which  may  be  the  case  of  Mars  (Bills  2006).  Fourth,  it  still  remains  a  matter  of 
controversy  as  to  whether  the  observed  pattern  of  small  craters  on  Mars  confirms  (Hartmann 

2007)  or  disproves  (Paige  et  al.  2007)  the  strong  climatic  variations  predicted  in  Ward  (1974, 
1979).  All  this  motivates  us  to  come  up  with  a  test  based  on  the  necessity  to  reconcile  the 
variations  of  spin  with  the  present  near-equatorial  positions  of  Phobos  and  Deimos.  The  fact 
that  both  moons  found  themselves  on  near-equatorial  orbits,  in  all  likelihood,  billions  of 
years  ago,1 *  and  that  both  are  currently  located  within  less  than  2  degrees  from  the  equator,  is 
surely  more  than  a  mere  coincidence.  An  elegant  but  sketchy  calculation  by  Goldreich  (1965) 
demonstrated  that  the  orbits  of  initially  near-equatorial  satellites  remain  close  to  the  equator 
of  date  for  as  long  as  some  simplifying  assumptions  remain  valid.  As  explained  in  Efroimsky 
(2004,  2005),  these  assumptions  are  valid  over  time  scales  not  exceeding  100  million  years, 
while  at  longer  times  a  more  careful  analysis  is  required.  Its  goal  will  be  to  explore  the  limits 
for  the  possible  secular  drift  of  the  satellite  orbits  away  from  the  evolving  equator  of  date. 
Through  comparison  of  these  limits  with  the  present  location  of  the  Martian  satellites,  we 
shall  be  able  to  impose  restrictions  upon  the  long-term  spin  variations  of  Mars.  If,  however,  it 
turns  out  that  near-equatorial  satellites  can,  in  the  face  of  large  equinoctial-precession  varia¬ 
tions,  remain  for  billions  of  years  close  to  the  moving  equator  of  date,  then  we  shall  admit  that 
Mars’  equator  could  indeed  have  precessed  through  billions  of  years  in  the  manner  predicted 
by  the  Colombo  approximation. 

The  second  motivation  for  our  study  comes  from  the  ongoing  discussion  of  whether  the 
Martian  satellites  might  have  been  captured  at  high  inclinations,  their  orbits  having  gradually 
approached  the  equator  afterwards.  While  a  comprehensive  check  of  this  hypothesis  will  need 
a  more  detailed  model — one  that  will  include  Mars’  triaxiality,  the  tidal  forces  (Lainey  et  al. 

2008) ,  and  perhaps  other  perturbations — the  first,  rough  sketch  of  this  test  can  be  carried  out 
with  only  J2 ,  the  Sun,  and  the  equinoctial  precession  taken  into  account.  We  perform  such  a 
rough  check  for  a  hypothetical  satellite  that  has  all  the  parameters  of  Deimos,  except  that  its 
initial  inclination  is  89°. 


1  Phobos  and  Deimos  give  every  appearance  of  being  captured  asteroids  of  the  carbonaceous  chondritic  type, 

with  cratered  surfaces  older  than  ~109  years  (Veverka  1977;  Pang  et  al.  1978;  Pollack  et  al.  1979;  Tolson  et  al. 

1978).  If  they  were  captured  by  gas  drag  (Burns  1972, 1978),  this  must  have  occurred  early  in  the  history  of  the 

solar  system  while  the  gas  disk  was  substantial  enough.  Kilgore,  Burns,  and  Pollack  (1978)  have  demonstrated 
numerically  that  a  gas  envelope  extending  to  about  ten  Martian  radii,  with  a  density  of  5  x  10-5  g/cm3 * * * * *  at  the 
Martian  surface,  could  have  been  capable  of  capturing  satellites  of  radii  about  10  km.  At  that  stage  of  planetary 
formation,  the  spin  of  the  forming  planet  would  be  perpendicular  to  the  planet’s  orbit  about  the  Sun — i.e.,  the 

obliquity  would  be  small  and  the  gas  disk  would  be  nearly  coplanar  with  the  planetary  orbit.  Energetically,  a 

capture  would  likely  be  equatorial.  This  is  most  easily  seen  in  the  context  of  the  restricted  three-body  problem. 
The  surfaces  of  zero  velocity  constrain  any  reasonable  capture  to  occur  from  directions  near  the  inner  and 

outer  collinear  Lagrange  points  (Szebehely  1967;  Murison  1988),  which  lie  in  the  equatorial  plane.  Also,  a 

somewhat  inclined  capture  would  quickly  be  equatorialised  by  the  gas  disk.  If  the  capture  inclination  is  too 
high,  the  orbital  energy  is  then  too  high  to  allow  a  long  enough  temporary  capture,  and  the  object  would 
hence  not  encounter  enough  drag  over  a  long  enough  time  to  effect  a  permanent  capture  (Murison  1988). 
Thus,  Phobos  and  Deimos  were  likely  (in  as  much  as  we  can  even  use  that  term)  to  have  been  captured  into 
near-equatorial  orbits. 
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1.3  Mathematical  tools 

The  first  steps  toward  the  analytical  theory  of  orbits  about  a  precessing  and  nutating  Earth  were 
undertaken  almost  half  a  century  ago  by  Brouwer  (1959),  Proskurin  and  Batrakov  (1960), 
and  Kozai  (1960).  The  problem  was  considered,  in  application  to  the  Martian  satellites,  by 
Goldreich  (1965)  and,  in  regard  to  a  circumlunar  orbiter,  by  Brumberg  et  al.  (1971).  The 
latter  two  publications  addressed  the  dynamics  as  seen  in  a  non-inertial  frame  co-precessing 
with  the  planet’s  equator  of  date.  The  analysis  was  carried  out  in  terms  of  the  so-called 
“contact”  Kepler  elements,  i.e.,  in  terms  of  the  Kepler  elements  satisfying  a  condition  differ¬ 
ent  from  that  of  osculation.  Modeling  of  perturbed  trajectories  by  sequences  of  instantaneous 
ellipses  (or  hyperbolae)  parameterised  with  such  elements  is  sometimes  very  convenient 
mathematically  (Efroimsky  2006c).  However,  the  physical  interpretation  of  such  solutions 
is  problematic,  because  instantaneous  conics  defined  by  nonosculating  elements  are  non¬ 
tangent  to  the  trajectory.  Though  over  restricted  time  scales  the  secular  parts  of  the  contact 
elements  may  well  approximate  the  secular  parts  of  their  osculating  counterparts  (Efroimsky 
2004,  2005),  the  cleavage  between  them  may  grow  at  longer  time  scales.  This  is  the  reason 
why  a  practically  applicable  treatment  of  the  problem  must  be  performed  in  the  language  of 
osculating  variables. 

1.4  The  plan 

The  analytical  theory  of  orbits  about  a  precessing  oblate  primary,  in  terms  of  the  Kepler  ele¬ 
ments  defined  in  a  co-precessing  (i.e.,  related  to  the  equator  of  date)  frame,  was  formulated  in 
Efroimsky  (2006a,  b)  where  the  planetary  equations  were  approximated  by  neglecting  some 
high-order  terms  and  averaging  the  others.  This  way,  from  the  exact  equations  for  osculating 
elements,  approximate  equations  for  their  secular  parts  were  obtained.  We  shall  borrow  those 
averaged  Lagrange-type  planetary  equations,  shall  incorporate  into  them  the  pull  of  the  Sun, 
and  shall  numerically  explore  their  solutions.  This  will  give  us  a  method  that  will  be  semi- 
analytical  and  seminumerical.  We  shall  then  apply  it  to  a  particular  setting — evolution  of  a 
Martian  satellite  and  its  reaction  to  the  long-term  variations  of  the  spin  state  of  Mars.  Our 
goal  will  be  to  explore  whether  the  spin- axis  variations  predicted  for  a  rigid  Mars  permit  its 
satellites  to  remain  close  to  the  equator  of  date  for  hundreds  of  millions  through  billions  of 
years.  In  case  the  answer  to  this  question  turns  out  to  be  negative,  it  will  compel  us  to  seek 
nonrigidity-caused  restrictions  upon  the  spin  variations.  Otherwise,  the  calculations  of  the 
rigid-Mars  inclination  variations  will  remain  in  force  (and  so  will  the  subsequent  calculations 
of  Mars  obliquity  variations);  this  way,  the  theory  of  Ward  (1973,  1974,  1979,  1982),  Touma 
and  Wisdom  (1993,  1994),  and  Laskar  and  Robutel  (1993)  will  get  a  model-independent 
confirmation. 


2  Semianalytical  treatment  of  the  problem 

To  understand  the  evolution  of  a  satellite  orbit  about  a  precessing  planet,  it  is  natural  to  model 
it  with  elements  defined  in  a  coordinate  system  associated  with  the  equator  of  date,  i.e.,  in 
a  frame  co-precessing  (but  not  co-rotating)  with  the  planet.  A  transition  from  an  inertial 
frame  to  the  co-precessing  one  is  a  perturbation  that  depends  not  only  upon  the  instantaneous 
position  but  also  upon  the  instantaneous  velocity  of  the  satellite.  It  has  been  demonstrated 
by  Efroimsky  and  Goldreich  (2004)  that  such  perturbations  enter  the  planetary  equations 
in  a  nontrivial  way:  not  only  do  they  alter  the  disturbing  function  (which  is  the  negative 
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Hamiltonian  perturbation)  but  also  they  endow  the  equations  with  several  extra  terms  that  are 
not  parts  of  the  disturbing  function.  Some  of  these  nontrivial  terms  are  linear  in  the  planet’s 
precession  rate  /x,  some  are  quadratic  in  it;  the  rest  are  linear  in  its  time  derivative  /x.  The 
inertial-forces-caused  addition  to  the  disturbing  function  (i.e.,  to  the  negative  Hamiltonian 
perturbation)  consists  of  a  term  linear  and  a  term  quadratic  in  /x.  (See  formulae  (53-54) 
in  Efroimsky  (2004,  2005)  or  formulae  (1)  and  (6)  in  Efroimsky  (2006a,  b).)  The  essence 
of  approximation  elaborated  in  Ibid,  was  to  neglect  the  quadratic  terms  and  to  substitute 
the  terms  linear  in  fi  and  /x  with  their  secular  parts  calculated  with  precision  up  to  e 3, 
inclusively. 

2. 1  Equations  for  the  secular  parts  of  osculating  elements  defined  in  a  co-precessing 
reference  frame 

We  shall  begin  with  five  Lagrange-type  planetary  equations  for  the  secular  parts  of  the  orbital 
elements  defined  in  a  frame  co-precessing  with  the  equator  of  date.  These  equations,  derived 
in  Efroimsky  (2004,  2005),  have  the  following  form: 

da  /,  2\V2 

—  =  -2 — a  (1  -  e  )  ,  (1) 

at  n 


de 

dt 


5  A_l  ,, 

- e  1 
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2 


n  h  / pe\2{5  2  •  1  \ 
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cos  i  I  .  /  3f  \  \ 
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in  i  \  \  dQjl 


(4) 


d£l  3  / Pe\2  COS  i  fin 

—  =  — n  .h  (  —  I  - w - 

dt  2  \a  /  (i  _  e2 y  sin  i 

+  2  12  1/2  ■  .  (a  (~f  X  (5) 

n  az{\  —  elyil  sin  i  \  \  di )  / 

The  number  of  equations  is  five,  because  one  element,  M0,  was  excluded  by  averaging 
of  the  Hamiltonian  perturbation  and  of  the  inertial  terms  emerging  in  the  right-hand  sides. 

In  the  equations,  n  =  Jg  (mprimary  +  m secondary )  A*3,  while  f  (V,  a,  e,  i,  co,  £2,  M0)  is  the 
implicit  function  that  expresses  the  unperturbed  two-body  dependence  of  the  position  upon 
the  time  and  Keplerian  elements.  Vector  fi  denotes  the  total  precession  rate  of  the  plane¬ 
tary  equator  (including  all  spin  variations — from  the  polar  wander  and  nutations  through  the 
Chandler  wobble  through  the  equinoctial  precession  through  the  longest-scale  spin  varia¬ 
tions  caused  by  the  other  planets’  pull),  while  /xi,  /X2,  /X3  stand  for  the  components  of  /x  in 
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a  co-precessing  coordinate  system  x,  y,  z,  the  axes  x  and  y  belonging  to  the  equator-of-date 
plane,  and  the  longitude  of  the  node,  £2,  being  reckoned  from  x: 

IL  =  /z ix  +  /X2y  +  /X3Z,  where  /zi  =  Ip,  fi 2  =  hp  sin  Ip,  /Z3  =  hp  cos  Ip,  (6) 

Ip,  hp  being  the  inclination  and  the  longitude  of  the  node  of  the  equator  of  date  relative  to 
that  of  epoch,  and  a  dot  standing  for  a  time  derivative.  The  quantity  /z_l  is  a  component  of  / 1 
directed  along  the  instantaneous  orbital  momentum  of  the  satellite,  i.e.,  perpendicular  to  the 
instantaneous  osculating  Keplerian  ellipse.  This  component  is  expressed  with 

/zj_  =  ft  •  w  =  sin  i  sin  £2  —  /Z2  sin  i  cos  £2  +  /z 3  cos  i,  (7) 


the  unit  vector 


w  =  x  sin  i  sin  £2  —  y  sin  i  cos  £2  +  z  cos  i  (8) 

standing  for  the  unit  normal  to  the  instantaneous  osculating  ellipse.  Be  mindful  that  /zj_  is 
defined  not  as  d(/JL  •  w )/dt  but  as 

=  [i  •  w  =  fl\  sin  i  sin  £2  —  1I2  sin  i  cos  £2  +  /Z3  cos  i.  (9) 

The  quantity  fin  is  a  component  pointing  within  the  satellite’s  orbital  plane,  in  a  direction 
orthogonal  to  the  line  of  nodes  of  the  satellite  orbit  relative  to  the  equator  of  date: 

fin  =  —fi\  sin  £2  cos  i  +  /z 2  cos  £2  cos  i  +  /Z3  sin  i 

=  —Ip  sin  £2  cos  i  +  hp  sin  Ip  cos  £2  cos  i  +  hp  cos  Ip  sin  i.  (10) 

Its  time  derivative  taken  in  the  frame  of  reference  co-precessing  with  the  equator  of  date  is: 

fin  =  —fi\  sin  £2  cos  i  +  {12  cos  £2  cos  i  +  /Z3  sin  i 

=  —Ip  sin  £2  cos  i  +  ( hp  sin  Ip  +  hpIp  cos  Ip )  cos  £2  cos  i 

+  (hp  cos  Ip  —  hplp  sin  Ip)  sin  i.  (11) 

As  shown  in  Efroimsky  (2006a,  b),  the  jx -dependent  terms,  emerging  in  Eqs.  1-5,  are 
expressed  with 

<z2 

—  {/zi  [—(2+3  e 2)  cos  ^2+5  e 2  (cos  ^2cos2<z>  —  sin  £2  sin2<x>cos  /)] 

+  A2  [—  (2  +  3  e 2)  sin  ^2  +  5  e 2  (sin  ^2  cos  2 co  +  cos  £2  sin  2 00  cos  /)] 
+  /Z3  [5^2  sin  2 00  sin  i  ]}  02) 


/t 


(2  +  3^2)  (/z  1  sin  i  sin  £2 


/Z2  sin  i  cos  £2  +  /Z3  cos  i), 


=  —  {/zi  sin  i  [—  (2  +  3  e 2)  sin  £2  cos  i 

+  5^2  (cos  £2  sin  loo  +  sin  £2  cos  2 00  cos  /)] 

+  /I2  sin  i  [  (2  +  3  e 2)  cos  ^2  cos  i 
+  5e2  (sin  ^2  sin  loo  —  cos  ^2  cos  2 00  cos  /)] 

—  /Z3  [(2  +  3  e 2)  (2  —  sin2  /)  +  5  e1  sin2  i  cos  2<z>]} 


(13) 


(14) 
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To  integrate  Eqs.  1-5,  with  expressions  (12-14)  inserted  therein,  we  shall  need  to  know,  at 
each  step  of  integration,  the  components  of  fi  and  (i. 

2.2  Calculation  of  the  components  of  /x  and  fi 

At  each  step  of  our  integration,  the  components  of  ji  and  (i  will  be  calculated  in  the  Colombo 
approximation.  Physically,  the  essence  of  this  approximation  is  two-fold:  first,  the  solar  torque 
acting  on  the  planet  is  replaced  by  its  average  over  the  year;  and,  second,  the  precessing  planet 
is  assumed  to  be  always  in  its  principal  spin  state.  While  a  detailed  development  (based  on 
the  work  by  Colombo  1966)  may  be  found  in  the  Appendix  to  Efroimsky  (2006a,  b),  here 
we  shall  provide  a  concise  list  of  resulting  formulae  to  be  used. 

The  components  of  fi  are  connected,  through  the  medium  of  (6),  with  the  inclination  and 
the  longitude  of  the  node  of  the  moving  planetary  equator,  I p  and  hp,  relative  to  some  equator 
of  epoch.  These  quantities  and  their  time  derivatives  are  connected  with  the  unit  vector  k 
aimed  in  the  direction  of  the  major-inertia  axis  of  the  planet: 

/v  T 

k  =  (sin Ip  sin/z^,  —  sin Ip  cos hp,  cos Ip)  .  (15) 

This  unit  vector  and  its  time  derivative 

dk 

—  =  [Ip  cos  Ip  sin hp  +  hp  sin  Ip  cos  hp, 

\T 

— Ip  cos  Ip  cos  hp  +  hp  sin  Ip  sin/z^,  —  Ipsinlp)  ,  (16) 

depend,  through  the  Colombo  equation 

f  =  “(*■*)  (fcxA)>  (17) 

upon  the  unit  normal  to  the  planetary  orbit, 

^  T 

n  =  (sin  Iorb  sin  Q0rb,  —  sin  Iorb  cos  £20rb,  cos  Iorb)  ,  (18) 

Slorb  and  Iorb  being  the  node  and  inclination  of  the  orbit  relative  to  some  fiducial  fixed  plane, 
and  a  being  a  parameter  proportional  to  the  oblateness  factor  fo.  In  our  computations,  we 
employed  the  present-day2  value  of  a — see  Table  5  below.  We  chose  this  value  because  it  was 
the  one  used  by  Ward  (1973,  1974),  and  we  wanted  to  make  sure  that  our  plot  for  obliquity 
evolution,  Fig.  3,  coincided  with  that  of  Ward  (1974). 

To  find  the  components  of  /jl ,  one  must  know  the  time  evolution  of  hp  and  Ip ,  which  can 
be  determined  by  solving  a  system  of  three  differential  equations  (17),  with  £20rb  and  Iorb 
being  some  known  functions  of  time.  These  functions  may  be  computed  via  the  auxiliary 
variables 


Q  —  Sin  I 0rb  Sin  ^2iorb ,  P  —  Sin  IQrb  COS  ^^orbi 


(19) 


2  An  accurate  treatment  shows  that  due  to  the  precession  of  the  Martian  orbit  a  exhibits  quasi-periodic 
variations  of  about  3%  over  long  time  scales.  Neglecting  this  detail  in  the  current  work,  we  shall  take  it  into 
account  at  the  further  stage  of  our  project  (Lainey  et  al.  2008). 
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Table  1  Numerical  values 
used  in  Ward’s  model  of  the 

j 

Ni 

s'j  [arcsec/yr] 

S'j  [deg] 

inclination  and  node  of  the 
Martian  orbit 

1 

0.0018011 

-5.201537 

272.06 

2 

0.0018012 

-6.570802 

210.06 

3 

-0.0358910 

-18.743586 

147.39 

4 

0.0502516 

-17.633305 

188.92 

5 

0.0096481 

-25.733549 

19.58 

6 

-0.0012561 

-2.902663 

207.48 

7 

-0.0012286 

-0.677522 

95.01 

whose  evolution  will  be  given  by  the  formulae: 


X  nj  sin  v 

7  =  1 

oo 

(20) 

X  nj  cos  ( 

7  =  1 

s'jt  +  8j)  . 

(21) 

The  choice  of  amplitudes,  frequencies,  and  phases  as  in  Table  1  will  make  equations  (19-21) 
render  Q0rb  and  Iorb  relative  to  the  solar  system’s  invariable  plane  (Ward  1974). 3 

The  development  by  Ward  (1974)  is  limited  in  precision.  A  more  accurate  treatment  was 
offered  by  Laskar  (1988).  At  the  future  stages  of  our  project,  when  developing  a  detailed 
physical  model  of  the  satellite  motion,  we  shall  employ  Laskar’ s  results.  In  the  current  paper, 
we  are  checking  if  the  orbital  averaging  of  the  precession-caused  terms  is  permissible  at  large 
time  scales.  Hence,  for  the  purpose  of  this  check,  as  well  for  the  qualitative  estimate  of  the 
long-term  behaviour  of  the  satellites,  we  need  a  realistic,  not  necessarily  highly  accurate, 
scenario  of  the  precession  variations. 

2.3  The  Goldreich  approximation 

The  above  semianalytical  treatment  not  only  yields  plots  of  the  time  dependence  of  the  mean 
elements  but  also  serves  as  a  launching  pad  for  analytical  approximations.  For  example,  an 
assumption  of  a  and  e  being  constant,  and  a  neglect  of  the  /(-dependent  terms  in  (4-5),  as 
well  as  of  the  term  /in/ sin  i  in  (5),  gives  birth  to  the  Goldreich  (1965)  approximation: 


-=0. 

dt 

^=0, 

dt 

di  . 

—  =  —fji  i  cos  £2  —  ii2  sin  £2, 
dt 


(22) 

(23) 

(24) 


3  Ward  (1974)  has  calculated  the  values  of  the  coefficients  N,  s' ,  8  based  on  the  work  of  Brouwer  and  van 
Woerkom  (1950)  who  had  calculated  the  values  of  p  and  q  relative  to  the  ecliptic  plane  of  1950.  Ward  (1974) 
transformed  Brouwer  and  van  Woerkom’s  values  to  the  solar  system’s  invariable  plane.  Be  mindful  that,  no 
matter  what  the  reference  plane,  both  Ward  (1974)  and  Brouwer  and  van  Woerkom  (1950)  used  the  same 
epoch,  J1950. 
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d£l 

dt 


cos  i 

(l-e2)2’ 


dco 

dt 


■*(?) 


Pe\2  5  cos2  i 


(1  -  £2)2 


1  cos  i 

H - : — : - 

sinz 


the  equinoctial  precession  being  assumed  uniform: 


iP  =  0 

hp  =  —a  cos  Ip 
Ml  =  0 

M2  =  hp  sin  Ip 
M3  =  hP  cos  Ip. 


(25) 

(26) 


(27) 

(28) 

(29) 

(30) 

(31) 


For  the  details  of  Goldreich’s  approximation  see  also  Subsect.  3.3  in  Efroimsky  (2005). 


2.4  The  gravitational  pull  of  the  Sun 


Reaction  of  a  satellite  on  the  planetary-equator  precession  is,  in  a  way,  an  indirect  reaction 
of  the  satellite  on  the  presence  of  the  Sun  and  the  other  planets.  Indeed,  the  pull  of  the  other 
planets  makes  the  orbit  of  our  planet  precess,  which  entails  variations  in  the  Sun-produced 
gravitational  torque  acting  on  the  planet.  These  variations  of  the  torque,  in  their  turn,  lead  to 
the  variable  equinoctial  precession  of  the  equator,  precession  “felt”  by  the  satellite.  It  would 
be  unphysical  to  consider  this,  indirect  effect  of  the  Sun  and  the  planets  upon  the  satellite, 
without  taking  into  account  their  direct  gravitational  pull.  In  this  subsection,  we  shall  take 
into  account  the  pull  of  the  Sun,  which  greatly  dominates  that  of  the  planets  other  than  the 
primary. 

In  what  follows,  m  and  m'  will  be  the  masses  of  the  satellite  and  the  Sun,  correspondingly, 
r  and  r'  will  stand  for  the  planetocentric  positions  of  the  satellite  and  the  Sun,  S  will  signify 
the  angle  between  these  vectors.  Then  the  Sun-caused  perturbing  potential  Rs,  acting  on  the 
satellite,  will  assume  the  form  of 


,/  1  r'-r\  Gm'  /  r'  rcosS\ 

Rs  Gm  (|r,  _r|  r/3  j  \|r'-r|  r'  ) 


(32) 


This  can  be  expanded  in  a  usual  manner  over  the  Legendre  polynomials  of  the  first  kind. 
Since  r'  r,  we  shall  take  only  the  first  term  in  the  series:4 


Rs 


Gm 

r 


ti'  T  r  /  r  \ 2 

—  1  +  —  Pi(cos£)  +  j  P2(cosS) 


r  cos  S' 


(33) 


As  the  term  Gm' /r'  is  not  dependent  of  the  satellite’s  elements,  one  has  only  to  consider 
the  restricted  potential 


i  Gm;  T/r  \  2  1  n'2a2  /a' \ 

=  —  |b)  —  (77) 


-  (3  cos2  S  —  1), 


(34) 


n'  and  a'  being  the  mean  motion  and  the  semi-major  axis  of  the  Sun. 

To  obtain  the  Lagrange-type  equations  for  the  third-body  perturbation,  we  must  first 
derive  an  expression  for  the  angle  S.  To  that  end,  define  the  directional  cosines  £  =  P  •  r'  and 


4  This  is  justified  since  the  next  term  in  the  Legendre  series  is  about  2  orders  of  magnitude  smaller  than  the 
potential  variations  generated  by  the  precession. 
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0  =  Q  •  r',  where  F  is  a  unit  vector  pointing  from  the  planet  toward  the  Sun,  while  P  and  Q 
are  unit  vectors  of  a  perifocal  coordinate  system  associated  with  the  osculating  orbital  plane 
of  the  satellite,  with  P  pointing  to  the  instantaneous  periapse  and  Q  being  orthogonal  to  P. 
Assuming  that  the  planet’s  orbit  about  the  Sun  is  circular,  we  arrive  at 

§  =  costucos(£2  —  M')  —  cos  i  sin  co  sin(£2  —  M'),  (35) 

0  =  —  sin  co  cos(£2  —  Mr)  —  cos  i  cos  co  sin(£2  —  M'),  (36) 

where  M'  is  the  mean  anomaly  of  the  Sun  in  the  planetocentric  frame.  With  aid  of  these 
relations,  the  angle  S  may  be  written  down  as 

cos  S  =  §  cos  v  +  6  sin  v,  (37) 


v  being  the  true  anomaly  of  the  satellite  in  the  planetocentric  coordinate  system.  Substituting 
(35)  and  (36)  into  (37),  and  averaging  over  the  satellite’s  mean  anomaly,  we  arrive  at  the 
Lagrange-type  equations: 


da 

dr 

de 

dr 


di 

dr 


dco 

dr 


d& 

dr 


(38a) 


10  e 


a/  1  —  e2  |^si 


sin2  i  sin  2co  +  (2  —  sin2  i)  sin  2 co  cos 


+  2  cos  i  cos  2 co  sin  2 £2 


] 


2  sin  i 


5e 2  cos  i  sin  2co(l  —  cos  2£2) 


Vl  —  e2 

-[2  +  e2(3  +  5  cos2&>)]  sin  2^| 

|4  +  e2  —  5  sin2  i  +  5 (sin2  i  —  e 2)  cos  2 co  +  5 (e2  —  2) 


(38b) 


(38c) 


Vl  —  e2 


i  sin  2 co  sin  2Q  +  [5(2  —  e2  —  sin2  i)  cos  2oo  —  2  —  3  e2  +  5  sin2  /] 


x  cos  2^2 


vT^ 


-5 e  sin  2 co  sin  2f2 


|  [2  +  e2(3  —  5  cos  2co)]  (1  —  cos  2^2)  cos  i 

i 


(38d) 


(38e) 


where  we  used  the  following  set  of  notations: 


r  =  /3n(t 


to ), 


3  m'a3 
16  Ma'3 


(39) 


16n  /  M\ 

3  n'  \  +  m' ) 


(40) 
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Q  =  Q-X', 


(41) 


A/  being  the  mean  longitude  of  the  Sun  in  the  planet’s  frame,  and  M  being  the  mass 
of  Mars. 

An  additional  averaging  can  be  performed  over  the  motion  of  the  planet  about  the  Sun. 
Mathematically,  this  is  the  same  as  averaging  over  the  motion  of  the  Sun  about  the  planet — in 
both  cases  averaging  over  A/  is  implied.  This  averaging  (Innanen  et  al.  1997)  will  simplify 
(38)  into 


da 

dt 

de 

dt 

di 

dt 

doo 

dt 

dQ 

dt 


0 


10ey/ 1  —  e2/3n  sin2  i  sin  2 o) 

\0e2  sin  i 

rfin  cos  i  sin  2 co 


VT^‘ 

2 

/3n  [4  +  e2  —  5  sin2  i  +  5 (sin2  i  —  e 2)  cos  2co] 


Vl  —  e2 


vT^ 


/3n  [2  +  £2(3  —  5  cos2&>)]  cos  i 


(42a) 

(42b) 

(42c) 

(42d) 

(42e) 


It  is  important  to  note  that  the  calculations  leading  to  equations  (38),  and  to  their  double- 
averaged  version,  (42),  were  performed  in  the  Martian-orbital,  and  not  Martian-equatorial 
frame,  without  taking  either  the  Martian  obliquity  or  precession  into  consideration.  Had  we 
taken  into  account  the  precession,  we  would  get,  on  the  right-hand  side  of  (42)  resonances 
between  the  motion  of  the  Sun  relative  to  the  planet  and  the  equinoctial  precession.  Since  the 
time  scale  of  the  former  exceeds,  by  orders  of  magnitude,  the  time  scale  of  the  latter,  we  may 
safely  omit  such  resonances.  This  justifies  our  neglect  of  the  frame  precession  in  the  above 
calculation. 

However,  the  omission  of  the  obliquity  may  have  a  serious  effect  on  the  results.  Thus,  we 
shall  generalize  equations  (38)  so  as  to  include  the  effect  of  the  solar  inclination  and  node 
in  a  Martian-centric  frame.  This  generalized  model  based  on  the  celebrated  works  by  Kozai 
(1959)  and  Cook  (1962)  gives  us: 


da 

dt 

de 

dt 

di 

dt 

do 

dt 


d£l 

dt 


0 


I5n'2e^/l  —  e 2 


4  n 


[2 AB  cos(2<z>)  —  (A2  —  B 2)  sin(2<z>)] 


3  n'2C 


4  n\[Y 


—(2  cos  /  + 


{A  [2  +  3^2  +  5^2  cos(2<z>)]  +  5 Be2  sin(2<z>)] 


3n/2Vl  —  e 2 

2  n 


x  5 AB  sin(2<z>)  +  -  (A2  —  B2)  cos(2<z>)  —  1  + 


3(A2  +  B2) " 


+ 


15 n'2a(A  cos  o)  +  B  sin  oo) 


4nea' 


1  -  ~(AZ  +  Bl) 


3n'2C 


4n  VT" 


—  [5Ae2  sin(2 co)  +  B  [2  +  3<?2  -  5e2  cos(2w)]} 


(43a) 

(43b) 

(43c) 


(43d) 

(43e) 
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where 


A  =  cos(£2  —  Qf)  cos(A/)  +  cos (/')  sin(A/)  sin(£2  —  Q')  (44a) 

B  =  cos  i  [—  sin(£2  —  £Y)  cos(A/)  +  cos(Y)  sin(A/)  cos(£2  —  Q')] 

+  sin  i  sin  (/')  sin  (A.7)  (44b) 

C  =  sin  i  [cos  A.7  sin(£2  —  £Y)  —  cos(Y)  sin  (A/)  cos(£2  —  £2')] 

+  cos  i  sin  i  ’  sin  X'  (44c) 


Here,  i'  and  £2'  are  the  inclination  and  right  ascension  of  the  ascending  node  of  the  Solar 
orbit  in  the  Martian  equatorial  frame,  respectively.  The  doubly  averaged  equations,  obtained 
after  averaging  over  the  Sun’s  mean  anomaly  for  a  single  period,  are  given  by 


da 

—  =  ° 
dt 

de 

dt 

di 


15 nae\J  1  —  e2  r  - 


4  n 


3  nL 


dt  4  n\J  l  —  e2 

dco  .  3  n/2V  1  —  e2 

—  =  —  Q  cos  i  4 - 

dt  In 


[2 AB  cos(2 co)  -  ( A 2  -  B2)  sin(2®)] 

[CA  [2  +  3e 2  +  5e2  cos(2<u)]  +  5 CBe2  sin(2<u)} 


_  5  _  _  3(A2  4-  B2) 

5 AB  sin(2 co)  +  -(A2  —  B2)  cos(2<u)  —  1  H - - - 


dQ 

dt 


3  n 


a 


(45a) 

(45b) 

(45c) 

(45d) 


An  V 1  —  e2  sin  i 


{5CA<?2  sin(2cu)  +  CB  [2  +  3^2  —  5e2  cos(2<u)]}  (45e) 


the  averaged  quantities  A2,  B2,  AB,  AC,  BC  being  given  by 

A2  =  [sj/Cjj,  -  0.5  s;,J  cQ  +  s;,sq/Sq/SqCq/Cq  +  0.5  -  0.5s(,cn,  (46a) 

B 2  =  {(s?,<4,  +  0.5s?/)  -  Sq/S2,S£2Cq'Cq  +  0.5s?, Cq,  -  0.5  +  c?,}  c? 

+  (Ci'CQCQ,  +  Cj/snsjy)  Sj/SjC,-  +  0.5s?,  (46b) 

AB  =  {s?,S£2'Cq'Cq  +  [  S?,  SqCq,  +  0.5s2,Sa]  Cq  -  0.5s2  Sft'CQ'}  C i 

+  0.5c ;/  (SqCq/  -  CqSq/)  S;/S;  (46c) 

AC  =  0.5  C ;/  (SqCq/  -  CqSq/)  Si/C i 

+  {-S?,Sq/Cq/Cq  +  [-s?,((sq)2Cq,)]  cq  +  0.5s-,sq/Cq/}  Si  (46d) 

BC  =  (Cj'CqCq/  +  Cj'SqSq/)  S;'C2 

+  {[4  (c£2')2  “  °-5  Cv)2]  CQ  -  S2Sq/SqCq/Cq  +  -0.5s?, Cq/  -  c2  +  0.5}  S/Ci 
-  0.5  (C//SqSq/  +  C;'CqCq/)  S//  (46e) 


where  we  have  used  the  compact  notation  cq  =  cos(-),  S(.)  =  sin(-).  To  calculate  the  trig¬ 
onometric  functions  of  &  and  i'  appearing  in  equations  (46),  we  shall  utilise  the  geometry 
rendered  Fig.  1.  The  figure  depicts  the  Martian  spin  axis,  k,  that  is  perpendicular  to  the  equa¬ 
tor  of  date,  and  the  normal  to  the  orbital  plane  (ecliptic  of  date),  n.  The  Martian  obliquity,  e, 
is  the  angle  between  these  two  vectors,  and  is  calculated  based  on  the  Colombo  formalism: 

cos  €  =  k  n  =  q  sx  +  p  sy  -\-  F  sz  (47) 
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Fig.  1  The  geometry  of  the  precessing  Mars.  The  Martian  spin  axis,  k,  is  perpendicular  to  the  equator  of  date, 
and  the  normal  to  the  orbital  plane  (ecliptic  of  date),  n.  The  Martian  obliquity,  e,  is  the  angle  between  these 
two  vectors 


where 

sx  =  sin  Ip  sin  hp,  sy  =  sin  Ip  cos  hp,  sz  =  cos  Ip,  F  =  f 1  —  p2  —  q2  (48) 
and  i'  =  €. 

The  angle  Qf,  lying  in  the  equator  of  date,  is  subtended  between  the  vectors  LON12  and 
LON23 .  The  first  of  these,  LON12,  points  along  the  line-of-nodes  obtained  by  the  intersection 
of  the  equator  of  date  and  the  invariable  plane.  This  vector  must  be  perpendicular  to  the  plane 
defined  by  k  and  z  (the  normal  to  the  invariable  plane),  so  that 

LON21  =  z  x  k  =  [sy,  sx,  of  (49) 

The  second  vector,  LON23,  is  aimed  along  the  line  of  nodes  rendered  by  the  intersection  of 
the  orbital  plane  and  the  equator  of  date.  Based  on  the  geometry  of  Fig.  1,  we  write: 

T 

LON23  =  k  x  n  =  [— syF  +  sz  p,  szq—sxF,  —sxp  +  syq]  (50) 


By  taking  the  scalar  product  of  LON23  and  LON21,  we  derive  the  direction  cosine: 

LON23  •  LON21  COS  Ip  COS  i'  —  COS  I  orb 


cos  £2'  = 


ILON23IILON21I 


sin  Ip  sin  /' 


(51) 


To  derive  an  expression  for  sin  £2',  we  shall  have  to  compute  an  auxiliary  vector,  j,  which  lies 
in  the  orbital  plane  and  is  normal  to  both  k  and  LON21  (cf.  Fig.  1): 

~  ry  T 

j  =  k  X  LON21  =  [-szsx,  szsy ,  sx+s  ] 


(52) 


The  direction  cosine  between  j  and  LON23  is  then 

(n  oA  LON23  •  j 

cos  ( - Q  )  = - — 

V2  )  |LON23 1  Ijl 

Upon  evaluating  Eq.  53  we  arrive  at 


(53) 
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P 


p  sin  hp 


(54) 


sin  V 


Finally,  the  calculation  of  the  orbital-elements’  evolution  is  performed  via  Lagrange-type 
planetary  equations,  whose  right-hand  sides  combine  those  of  (l)-(5)  and  (45). 

2.5  The  higher-order  harmonics,  the  gravitational  pull  of  the  planets, 
the  Yarkovsky  effect,  and  the  bodily  tides 

In  the  current  paper,  we  pursue  a  limited  goal  of  taking  into  account  the  oblateness  of  the 
planet,  its  nonuniform  equinoctial  precession,  and  the  gravitational  pull  of  the  Sun.  These 
three  items  certainly  do  not  exhaust  the  list  of  factors  influencing  the  orbit  evolution  of  a 
satellite. 

Among  the  factors  that  we  intend  to  include  into  the  model  at  the  further  stage  of  its  devel¬ 
opment  are  the  high-order  zonal  (^3,  J4,  J%)  and  tesseral  (C22)  harmonics  of  the  planet’s 
gravity  field,  as  well  as  the  gravitational  pull  of  the  other  planets — factors  whose  role  was 
comprehensively  discussed,  for  example,  by  Waz  (2004).  We  also  intend  to  include  the  bodily 
tides  (Efroimsky  and  Lainey  2007)  and  the  Yarkovsky  effect  (Nesvorny  and  Vokrouhlicky 
2007) — factors  whose  importance  increases  at  long  time  scales. 

3  Comparison  of  a  purely  numerical  and  a  semianalytical  treatment 
of  the  problem 

One  of  our  goals  is  to  check  the  applicability  limits  (both  in  terms  of  the  initial  conditions  and 
the  permissible  time  scales)  of  our  semianalytical  model  written  for  the  osculating  elements 
introduced  in  a  frame  co-precessing  with  the  equator  of  date.  This  check  will  be  performed  by 
an  independent,  purely  numerical,  computation  that  will  be  free  from  whatever  simplifying 
assumptions  (all  terms  kept,  no  averaging  performed.)  The  straightforward  simulation  will  be 
carried  out  in  terms  of  Cartesian  coordinates  and  velocities  defined  in  an  inertial  frame  of  ref¬ 
erence.  Both  the  semianalytical  calculation  of  the  elements  in  a  co-precessing  frame  and  the 
straightforward  numerical  integration  in  inertial  Cartesian  axes  will  be  carried  out  for  Deimos. 

3.1  Integration  by  a  purely  numerical  approach 

The  numerical  integration  of  Deimos’  orbit  can  be  performed  using  Cartesian  coordinates 
defined  relatively  to  the  Solar  system  invariable  plane.  As  we  also  have  to  compute  the 
Martian  polar  axis  motion,  there  are  two  vector  differential  equations  to  integrate  simulta¬ 
neously.  One  is  the  Newton  gravity  law: 


where  V  U  1 
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Here  0  and  r  =  (x,  y,  z)  denote,  correspondingly,  the  latitude  of  Deimos  relative  to  the 
Martian  equator  and  the  position  vector  of  Deimos  related  to  the  Martian  center  of  mass;  pe 
is  the  Martian  equatorial  radius;  M  and  m  stand  for  the  masses  of  Mars  and  Deimos,  respec¬ 
tively.  Angles  hp  and  Ip  are  the  longitude  of  the  node  and  the  inclination  of  the  planet’s 
equator  of  date  relative  to  the  invariable  plane  (see  Sect.  2.2).  Integration  in  this,  inertial, 
frame  offers  the  obvious  advantage  of  nullifying  the  inertial  forces. 

Table  2  gives  the  initial  conditions  for  our  simulation,  expressed  in  terms  of  the  Keplerian 
orbital  elements.  Table  3  presents  these  initial  conditions  in  a  more  practical  form,  i.e.,  in 
terms  of  the  Cartesian  positions  and  velocities  corresponding  to  the  said  elements.  A  tran¬ 
sition  from  the  Keplerian  elements  to  these  Cartesian  positions  and  velocities  is  a  two-step 
process.  First,  we  take  orbital  elements  defined  in  a  frame  associated  with  the  Martian  equa¬ 
tor  of  date  (i.e.,  in  a  frame  co-precessing  but  not  co-rotating  with  the  planet)  and  transform 
them  into  Cartesian  coordinates  and  velocities  defined  in  that  same  frame.  Then,  by  two 
successive  rotations  of  angles  —Ip  and  —hp,  we  transform  them  into  Cartesian  coordinates 
and  velocities  related  to  the  invariable  plane.  These  initial  positions  and  velocities  were  used 
to  begin  the  integration. 

At  each  step  of  integration  of  (55),  the  same  two  rotations  are  performed  on  the  compo¬ 
nents  VI/  given  by  (56).  As  mentioned  above,  to  afford  the  absence  of  inertial  forces  on  the 
right-hand  side  of  (55)  one  must  write  down  and  integrate  (55)  in  the  inertial  frame.  Since  the 
analytical  expressions  (56)  for  VI/  contain  the  latitude  0,  they  are  valid  in  the  co-precessing 
coordinate  system  and,  therefore,  need  to  be  transformed  to  the  inertial  frame  at  each  step.  To 
carry  out  the  transformation,  one  needs  to  know,  at  each  step,  the  relative  orientation  of  the 
Martian  polar  axis  and  the  inertial  coordinate  system.  The  orientation  is  given  by  the  afore 
mentioned  Colombo  model.  This  is  how  our  second  equation,  the  one  of  Colombo,  comes 
into  play: 

dk  ~  . 

—  =  a(k  •  n)(k  x  n).  (57) 

dt 

All  in  all,  we  have  to  integrate  the  system  (55-57).  Table  4  gives  the  initial  conditions 
used  for  integrating  (57),  while  Table  5  gives  the  numerical  values  used  for  the  parameters 


Table  2  The  orbital  elements 
values  taken  as  initial  conditions 

Parameters 

Numerical  values 

for  our  simulations 

a 

23459  km 

e 

0.0005 

i 

0.5  deg  and  89  deg 

Q 

10  deg 

CO 

5  deg 

M 

Odeg 

Table  3  The  initial  positions  and  velocities  used  for  Deimos.  The  first  two  rows  correspond  to  the  low- 
inclination  case  (0.5  degree);  the  last  two  rows  correspond  to  the  high-inclination  case  (89  degrees) 

Satellite 

X 

y 

z 

Position  km  (i  =  0.5) 
Velocity  km/s  (i  =  0.5) 
Position  km  (i  =  89) 
Velocity  km/s  (i  =  89) 

22648.3376439 

-0.349882011871 

22996.9921622 

-0.120115009144 

6068.52353055 

1.30576017694 
4091.20549954 
2.686751629968  xlO-3 

17.8332361962 
1.175229063323  xlO-2 
2043.25303109 

1.34652528539 
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Table  4  Initial  conditions  used 
for  the  integration  of  Eq.  57. 
These  values  were  calculated 
based  on  Ward  (1974) 


Table  5  Parameter  values 
used  in  our  simulations 


Parameters 

Numerical  values 

Ip(t0) 

25.25797549  deg 

hp{t o) 

332.6841708  deg 

Parameters 

Numerical  values 

Martian  mass  (GM) 

42830  km3  s-2 

h 

1960.45  x  10“6 

Equatorial  radius 

3397  km 

Deimos  mass 

0.091  x  10-3  km3  s-2 

Oi 

3.9735  x  10“5 *  rad/yr 

involved.  It  is  worth  noting  that  the  values  of  Table  4  were  calculated  based  on  Ward 
(1974). 

The  software  used  for  numerical  integration  of  the  system  (55-57)  is  called  NOE  (Numeri¬ 
cal  Orbit  and  Ephemerides),  and  is  largely  based  on  the  ideas  and  methods  developed  in  Lainey 
et  al.  (2004).  This  numerical  tool  was  created  at  the  Royal  Observatory  of  Belgium  mainly  for 
computations  of  the  natural  satellite  ephemerides.  It  is  an  A-body  code,  which  incorporates 
highly  sensitive  modelling  and  can  generate  partial  derivatives.  The  latter  are  needed  when 
one  wants  to  fit  the  initial  positions,  velocities,  and  other  parameters  to  the  observation  data. 
To  save  the  computer  time,  an  optimised  force  subroutine  was  built  into  the  code,  specifically 
for  integrating  the  above  equations.  This  appliance,  based  on  the  RA15  integrator  offered  by 
Everhart  (1985),  was  chosen  for  its  speed  and  accuracy.  During  the  integration,  a  variable 
step  size  with  an  initial  value  of  0.04  day  was  used.  To  control  the  numerical  error,  back  and 
forth  integrations  were  performed.  In  particular,  we  carried  out  a  trial  simulation  consisting 
of  a  thousand-year  forward  and  a  subsequent  thousand-year  back  integration.  The  satellite 
displacement  due  to  the  error  accumulated  through  this  trial  was  constrained  to  150  m.  Most 
of  this  150  m  difference  comes  from  a  numerical  drift  of  the  longitude,  while  the  numerical 
errors  in  the  computation  of  the  semi-major  axis,  the  eccentricity,  and  the  inclination  were 
much  lower.  These  errors  were  reduced  for  this  trial  simulation  to  only  10-5  km,  10-10, 
and  10“ 10  degree,  respectively.  This  provided  us  with  a  high  confidence  in  our  subsequent 
numerical  results. 

As  a  complement  to  the  said  back-and-forth  check,  the  energy-conservation  criterium  was 
used  to  deduce,  in  the  first  approximation,  an  optimal  initial  step-size  value  and  to  figure  out 
the  numerical  error  proliferation.  (It  is  for  this  energy-conservation  test  that  we  introduced 
a  non-zero  mass  for  Deimos.  Its  value  was  taken  from  Smith  et  al.  (1995).)  Applicability  of 
this  criterium  is  justified  by  the  fact  that  the  numerical  errors  are  induced  mostly  by  the  fast 
orbital  motion  of  the  satellite.5 


5  Although  the  planet-satellite  system  is  subject  to  an  external  influence  (the  solar  torque  acting  on  the  planet), 

over  short  time  scales  this  system  can  be  assumed  closed.  In  order  to  check  the  integrator  efficiency  and  to 

determine  an  optimal  initial  step  size,  we  carried  out  auxiliary  integrations  of  (55)-(56),  with  the  Colombo 

equation  (57)  neglected  and  with  the  energy  presumed  to  conserve.  These  several-thousand-year-long  trial  inte¬ 

grations,  with  the  energy-conservation  criterium  applied,  led  us  to  the  conclusion  that  our  integrator  remained 
steady  over  long  time  scales  and  that  the  initial  step  of  0.04  day  was  optimal.  Then  this  initial  step  size  was 
employed  in  our  integration  of  the  full  system  (55)-(57). 
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3.2  Integration  by  the  semianalytical  approach 

The  theory  of  satellite-orbit  evolution,  based  on  the  planetary  equations  whose  right-hand 
sides  combine  those  of  (1-5)  and  (45),  is  semianalytical.  This  means  that  these  equations  for 
the  elements’  secular  parts  are  derived  analytically,  but  their  integration  is  to  be  performed 
numerically.  This  integration  was  carried  out  using  an  %th -order  Runge-Kutta  scheme  with 
relative  and  absolute  tolerances  of  10-12.  Kilograms,  years,  and  kilometers  were  taken  as 
the  mass,  time,  and  length  units,  correspondingly. 


3. 2. 1  Technicalities 


To  integrate  the  planetary  equations,  one  should  know,  at  each  time  step,  the  values  of  Ip  and 
hp,  which  are  the  time  derivatives  of  the  inclination  and  of  the  longitude  of  the  node  of  the 
equator  of  date  with  respect  to  the  equator  of  epoch.  These  derivatives  will  be  rendered  by 
the  Colombo  equation  (17),  after  formulae  (15-16)  get  inserted  therein: 

Ip  =  —a  ^ q 2  sin  Ip  sin  hp  cos  hp  —  qp  sin  Ip  +2 pq  sin  Ip  cos2  hp 
—  p1  sin  Ip  cos  hp  sin  hp  +  qj  1  —  q2  —  p2  cos  Ip  cos  hp 

7 


and 


hp  =  —a 


—  pyj  1  —  q2  —  p2  cos  Ip  sin  hp 

[—q  +  2 q  cos2  Ip)  cos2  hp  —  2 q  cos2  Ip  +  q 


(58) 


(/?  —  2/7  cos2  Ip)  cos  hp  + 


sin  hr 


.  yi  -p2  -< 

sin  In 


+  (#2  —  P 2)  cos  h  cos2  hp  +  (— P2  ~  2 q2  +  l)  cos  Ip 


+ 


2qp  cos  Ip  cos3  hp  —2 pq  cos  hp  cos  Ip 
sin  hr, 


(59) 


Equations  (58)  and  (59)  are  then  integrated  (with  the  initial  conditions  Ip  (to)  and  hp  (to)  bor¬ 
rowed  from  Table  4)  simultaneously  with  the  planetary  equations  (1-5).  Through  formulae 
(6),  the  above  expressions  for  Ip  and  hp  yield  the  expressions  for  the  components  of  ii.  As 
can  be  seen  from  (12-14),  integration  of  the  planetary  equations  also  requires  the  knowledge 
of  the  derivatives  jl\  ,  /X2,  and  /X3  at  each  integration  step.  These  can  be  readily  obtained 
by  differentiating  (6).  The  resulting  closed-formed  expressions  for  (i\,  fi2,  M3  are  listed  in 
Appendix  A.  The  final  step  required  for  numerical  integration  of  the  planetary  equations  is 
substitution  of  formulae  (9-11),  with  the  initial  conditions  from  Table  2. 


3.2.2  The  plots  and  their  interpretation 

Figure  2  depicts  the  history  of  the  planetary  equator,  in  the  Colombo  approximation,  over  1 
Byr.  The  inclination  exhibits  long -periodic  oscillations  bounded  within  the  range  of  20.3  deg 
<  Ip  <  30.3  deg,  while  the  node  regresses  at  a  rate  of  hp  =  0.00202  deg  /yr.  The  obliquity, 
c,  varies  in  the  range  15.2  deg  <  c  <  35.5  deg.  A  magnified  view  of  the  obliquity  for  5  Myr 
is  shown  in  Fig.  3.  The  time  history  of  the  obliquity  closely  matches  the  results  reported  by 
Ward  (1974). 


Springer 


278 


P.  Gurfil  et  al. 


x  106 

5  - 1 - T - T - 1 - 

o  - — . : . ; . - 

-5  - 1 - 1 1 - 1 - 

0  2  4  6  8  10 


Time  [yr] 


x  10 


8 


Fig.  2  Evolution  of  the  Martian  inclination,  the  longitude  of  the  node  of  the  equator  of  date  relative  to  that  of 
epoch,  and  the  obliquity  over  1  Byr,  obtained  in  the  Colombo  approximation 


Time  [yr]  x1()6 

Fig.  3  Evolution  of  the  the  Martian  obliquity  for  5  Myr,  calculated  through  formula  (47).  This  curve  closely 
matches  the  result  of  Ward  (1974) 


Integration  of  our  semianalytical  model  gives  plots  depicted  in  Figs.  4  and  5,  for  a  low 
initial  inclination,  and  in  Figs.  6-7,  for  a  high  initial  inclination.  From  Fig.  4  we  see  that 
the  variable  equinoctial  precession  does  not  inflict  considerable  changes  upon  the  satellite’s 
inclination  relative  to  the  precessing  equator  of  date.  The  orbit  inclination  remains  bounded 
within  the  region  0.3  deg  <  i  <  2.5  deg.  This  means  that  the  “Goldreich  lock”  (inclination 
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Fig.  4  Evolution  of  the  inclination  (initially  set  to  0.5  degree)  and  of  the  longitude  of  the  node  of  Deimos 
over  1  Byr.  The  plot,  obtained  by  integration  of  the  semianalytical  model,  exhibits  inclination  locking  and  a 
uniform  regression  of  the  node 


“locking,”  predicted  by  the  Goldreich  (1965)  model  for  small  inclinations  and  for  uniform 
equinoctial  precession)  works  also  for  nonuniform  Colombo  precession  of  the  equator. 

However,  this  observation  no  longer  holds  for  large  initial  inclinations.  In  the  z'o  =  89  deg 
case,  it  is  clearly  seen  that  the  node  is  greatly  affected  by  the  presence  of  the  equinoctial  pre¬ 
cession.  The  equinoctial  precession  also  affects  the  magnitude  of  the  inclination  variations. 
This  case  exhibits  chaotic  dynamics  that  are  sensitive  to  any  additional  perturbing  inputs. 
Figure  6  shows  that  without  precession  the  inclination  of  Deimos’  orbit  varies  within  the 
range  of  84.5  deg  <  i  <  95  deg.  With  the  precession  included,  the  inclination  gains  about 
one  degree  in  amplitude,  varying  in  the  range  83.5  deg  <  i  <  96  deg.  This  effect  is  accentu¬ 
ated  when  examining  a  magnification  of  the  inclination  over  a  5  Myr  span,  as  shown  in  Fig.  6. 
The  chaotic  nature  of  the  inclination  is  clearly  seen.  The  irregular  dynamics  is  characterised 
by  chaotic  switches  between  the  maximum  and  minimum  inclination  values,  a  phenomenon 
referred  to  as  “crankshaft”,  to  be  further  discussed  in  the  sequel. 

Both  in  the  near-equatorial  case  (as  in  Fig.  5)  and  the  near-polar  case  (as  in  Fig.  7), 
variations  of  the  semimajor  axis  are  of  order  10-6%.  This  smallness  is  in  compliance  with 
formula  (44)  in  Efroimsky  (2006a,  b),  according  to  which  the  changes  of  a  generated  by 
the  variations  of  the  equinoctial  precession  rate  are  extremely  small.  (Formula  (38a)  from 
Subsect.  2.4  above  tells  us  that  the  direct  pull  of  the  Sun  exerts  no  influence  upon  a  at  all.) 

Similarly,  in  both  cases  (Figs.  5  and  7)  the  variations  of  eccentricity,  remain  small,  about 
10-3%.  While  formula  (44)  in  Efroimsky  (2006a,  b)  promises  to  e  only  tiny  variations  due  to 
precession,  our  formula  (38b)  from  Subsect.  2.4  above  gives  to  e  a  slightly  higher  variation 
rate,  rate  that  still  remains  insufficient  to  raise  the  quasiperiodic  changes  in  e  above  a  fraction 
of  percent. 
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Fig.  5  Evolution  of  the  semimajor  axis,  eccentricity,  and  argument  of  periapsis  of  Deimos  over  1  Byr.  (The 
inclination  was  initially  set  to  0.5  degree.)  Both  the  semimajor  axis  and  eccentricity  exhibit  quasiperiodic 
motion  about  their  initial  values.  (The  variations  of  the  semimajor  axis  and  eccentricity  are  so  small  that  it 
is  more  convenient  to  plot  a  —  a$  and  e  —  eg.)  The  plots  were  obtained  by  integration  of  the  semianalytical 
model 


The  plots  in  Figs.  5  and  7  depict  also  the  time  evolution  of  co.  The  line  of  apsides  steadily 
regresses  in  the  near-polar  case  and  steadily  advances  in  the  near-equatorial  case. 

To  examine  the  precision  of  our  semianalytical  model,  we  have  compared  the  results  of 
its  integration  with  the  results  stemming  from  a  purely  numerical  simulation  performed  in 
terms  of  inertial  Cartesian  coordinates  and  velocities  (see  Subsect.  3.1).  The  necessity  for  this 
check  was  dictated,  mainly,  by  the  fact  that  within  the  semianalytical  model  the  short-period 
terms  are  averaged  out,6  while  the  straightforward  numerical  integration  of  (55-56)  neglects 
nothing.  We  carried  out  comparison  of  the  two  methods  over  10  Myr  only.  The  outcomes, 
both  for  io  =  0.5  degrees  and  io  =  89  deg,  were  in  a  good  agreement.  As  an  example, 
the  top  plot  in  Fig.  8  shows  the  comparison  of  the  inclination  evolution  calculated  by  the 
semianalytical  and  purely  numerical  methods  over  10  Myr  in  the  case  of  io  =  0.5  deg.  The 
bottom  plot  in  Fig.  8  depicts  a  similar  comparison  for  the  case  of  io  =  89  deg.  Since  the 
inclination  exhibits  chaotic  behaviour,  comparing  the  semianalytical  and  purely  numerical 
calculations  point- by-point  (i.  e.,  computing  the  differences  between  these  two  signals)  would 
not  be  useful.  Therefore,  we  chose  to  compare  the  mean  and  standard  deviation  (STD)  of 
the  inclination,  in  the  semianalytical  and  purely  numerical  simulations.  We  also  compared 
the  extremum  values.  The  results  of  this  comparison  are  summarised  by  Table  6.  The  table 
quantifies  the  agreement  between  the  models  depicted  by  Fig.  8.  It  clearly  demonstrates  that 
the  two  simulations  agree  up  to  fractions  of  a  percent. 


6  We  remind  that  in  equations  (1-5)  the  exact  //.-dependent  terms  are  substituted  with  their  orbital  averages 
(12-14). 
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Fig.  6  Evolution  of  the  inclination  and  of  the  longitude  of  the  node  of  Deimos  over  1  Byr.  (The  inclina¬ 
tion  was  initially  set  to  89  degrees.)  The  plot,  obtained  by  integration  of  the  semianalytical  model,  exhibits 
inclination  variation  in  the  range  ±5  deg  and  a  chaotic  evolution  of  the  node;  if  precession  is  neglected,  the 
inclination  oscillations  are  smaller  in  magnitude.  The  chaotic  nature  of  the  inclination  variation  is  referred  to 
as  “crankshaft  chaos” 


All  in  all,  the  outcome  of  our  computations  is  two-fold.  First,  we  have  made  sure  that  the 
semianalytical  model  perfectly  describes  the  dynamics  over  time  scales  of,  at  least,  dozens 
of  millions  of  years.  Stated  differently,  the  short-period  terms  and  the  terms  of  order  0(fi2) 
play  no  role  over  these  time  spans.  Second,  we  have  made  sure  that  the  “Goldreich  lock” 
initially  derived  for  very  low  inclinations  and  for  uniform  equinoctial  precession,  works  well 
also  for  variable  precession,  though  for  low  initial  inclinations  only. 
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Fig.  7  Evolution  of  the  semi-major  axis,  eccentricity,  and  argument  of  periapsis  of  Deimos  over  1  Byr.  (The 
inclination  was  initially  set  to  89  degrees.)  The  semi-major  axis  exhibits  chaotic  behavior  and  the  eccentricity 
exhibits  long-periodic  motion  about  the  initial  value.  (The  variations  of  the  semimajor  axis  are  so  small  that  it 
is  more  convenient  to  plot  a  —  a0  rather  than  a.)  The  plots  were  obtained  by  integration  of  the  semianalytical 
model 


Table  6  Statistical  properties  of  the  inclination.  Comparison  of  the  results  of  the  semianalytical  and  the  purely 
numerical  computation 


z'o  [deg] 

Model 

STD  [deg] 

Mean  [deg] 

Max.  value  [deg] 

Min.  value  [deg] 

0.5 

Semianalytical 

0.60 

1.519 

2.45 

0.3063 

Cartesian 

0.601 

1.53 

2.465 

0.3056 

89 

Semianalytical 

3.10 

90.085 

95.9713 

84.027 

Cartesian 

3.09 

89.92 

95.9769 

84.0054 

3.2.3  Looking  for  trouble 

A  natural  question  arises  as  to  whether  the  considered  examples  are  representative.  One  may 
enquire  if,  perhaps,  there  still  exists  a  combination  of  the  initial  conditions  yielding  noticeable 
variations  of  the  satellite  orbit  inclination  during  the  primary’s  variable  equinoctial  preces¬ 
sion  over  vast  spans  of  time.  To  answer  this  question,  we  should  scan  through  all  the  possible 
combinations  of  initial  conditions,  to  identify  a  particular  combination  that  would  entail  a 
maximal  inclination  excursion  relative  to  the  initial  inclination  (Gurfil  et  al.  2002).  Stated 
more  formally,  we  should  seek  a  set  of  initial  conditions  {/z*0,  7*0,  z'q,  £2q,  cUq}  maximising 
the  objective  function  |  i(t)  —  z‘o|: 

{h*p0,  /*o,  «5,  £2o>  =  argmax|i(0  -  i0\  (60) 

rhp0,lp0 

io,Q,o,ojo 
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Fig.  8  Comparison  of  the  semi-analytical  model  to  a  purely  numerical  integration.  The  plot  shows  the  incli¬ 
nation  as  predicted  by  the  two  models.  The  top  plot  shows  the  case  of  the  initial  inclination  iq  =  0.5  deg.  In 
this  case,  the  standard  deviation  of  the  inclination  in  the  semianalytical  and  inertial  models  differs  by  0.175%, 
and  the  mean  value  by  0.77%.  The  bottom  plot  shows  the  case  of  the  initial  inclination  iq  =  89  deg.  In 
this  case,  the  differences  between  the  semianalytical  and  the  numerical  model  amount  to  0.4%  in  standard 
deviation  and  0.18%  in  the  mean  value 


This  problem  belongs  to  the  realm  of  optimisation  theory.  The  optimisation  space  is  consti¬ 
tuted  by  the  entire  multitude  of  the  permissible  initial  conditions.  The  sought  after  combina¬ 
tion  of  initial  conditions  will  be  called  the  optimal  set. 

In  this  situation,  the  traditional  optimisation  schemes  (such  as  the  gradient  search  or  the 
simplex  method)  may  fail  due  to  the  rich  dynamical  structure  of  our  problem — these  methods 
may  lead  us  to  a  local  extremum  only.  Thus,  the  search  needs  to  be  global.  It  can  be  carried 
out  by  means  of  Genetic  Algorithms  (GAs)  (Goldberg  1989;  Gurfil  et  al.  2002;  Gurfil  and 
Kasdin  2002a,  b).  These  are  wont  to  supersede  the  traditional  optimisation  procedures  in  the 
following  aspects.  First,  instead  of  directly  dealing  with  the  parameters,  the  GAs  employ 
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Table  7  Parameter  values  used 
for  the  GA  optimization 


Parameters 


Numerical  values 


Population  size 


30 
100 
16  bit 
0.99 
0.02 


Number  of  generations 
String  length 


Probability  of  crossover 
Probability  of  mutation 


codings  (usually,  binary)  of  the  parameter  set  (“strings,”  in  the  GA  terminology).  Second, 
instead  of  addressing  a  single  point  of  the  optimisation  space,  the  GAs  perform  a  search  inside 
a  population  of  the  initial  conditions.  Third,  instead  of  processing  derivatives  or  whatever 
other  auxiliary  information,  the  GAs  use  only  objective-function  evaluations  (“fitness  eval¬ 
uations”).  Fourth,  instead  of  deterministic  rules  to  reiterate,  the  GAs  rely  upon  probabilistic 
transition  rules.  Additional  details  on  the  particular  GA  mechanism  used  herein  can  be  found 
in  Appendix  B. 

A  GA  optimisation  was  implemented  using  the  parameter  values  given  in  Table  7. 

The  search  for  the  inclination-maximising  initial  conditions  resulted  in  the  following  set: 


/*0  =  72.5  deg,  h*p0  =  211.324  deg,  =  100.543  deg, 
=  111.538  deg,  o%  =  234.913  deg. 


(61) 


Thus,  the  initial  orbit  is  retrograde  and,  not  surprisingly,  near-polar.  The  resulting  time  his¬ 
tories  for  a  0.2  Byr  integration  are  depicted  in  Fig.  9,  for  i  and  £2.  In  both  cases  the  inclination 
amplitude  is  relatively  large:  The  inclination  varies  within  the  range  of  79  deg  <  i  <  102  deg. 

This  example  clearly  shows  that  the  equinoctial  precession  is  an  important  effect  for 
evolution  of  satellite  orbits.  As  shown  in  the  upper  pane  of  Fig.  9,  had  we  neglected  the 
precession,  the  magnitude  of  the  oscillations  would  be  about  twice  smaller.  The  inclusion 
of  the  precession  in  the  model  qualitatively  modifies  the  behavior,  inducing  large-magnitude 
chaotic  variations  of  the  inclination,  a  phenomenon  that  cannot  be  detected  without  including 
the  precession  alongside  the  oblateness  and  solar  gravity. 

4  Comparison  of  the  semianalytical  results  with  those  rendered 
by  Goldreich’s  model 

The  final  step  in  our  study  will  be  to  compare  the  semianalytical  model  to  Goldreich’s  approx¬ 
imation  (22-26).  To  that  end,  we  integrate  our  semianalytical  model  for  20  Myr,  using  the 
initial  conditions  from  Table  2  with  z'o  =  89  deg;  and  compare  the  outcome  with  that  resulting 
from  Goldreich’s  approximation  simulated  with  the  same  initial  conditions.  The  results  of  this 
comparison  are  depicted  in  Figs.  10  and  1 1 .  Specifically,  Fig.  10  compares  the  time  histories  of 
Ip  and  hp .  There  are  noticeable  differences  in  the  dynamics  of  Ip .  While  Goldreich’s  approx¬ 
imation  assumes  a  constant  Ip,  the  semianalytical  model  is  based  on  the  Colombo  calculation 
of  the  equinoctial  precession,  calculation  that  predicts  considerable  oscillations  within  the 
range  21  deg  <  Ip  <  30  deg.  Beside  this,  in  our  semianalytical  model  we  take  into  account 
the  direct  gravitational  pull  exerted  by  the  Sun  on  the  satellite.  All  this  entails  differences 
between  the  dynamics  predicted  by  our  semianalytical  model  and  the  dynamics  stemming 
from  the  Goldreich  approximation.  These  differences,  for  i  and  £2,  are  depicted  in  Fig.  11. 
In  Goldreich’s  model,  i  stays  very  close  to  the  initial  value:  88.27  deg  <  i  <  89.01  deg, 
a  behaviour  that  makes  the  essence  of  the  Goldreich  lock.  However,  in  the  more  accurate, 
semianalytical  model  we  have  84  deg  <  i  <  96  deg.  The  time  history  of  £2,  too,  reveals  that 
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Fig.  9  Evolution  of  the  inclination  and  longitude  of  the  node,  for  the  initial  conditions  that  entail  maximal 
variations  of  the  orbit  inclination.  The  plot,  obtained  by  integration  of  the  semianalytical  model,  demonstrates 
that  the  inclination,  initially  set  at  100.5  deg,  plunges  to  less  than  80  deg  and  then  returns  to  its  initial  value. 
The  switches  between  maximum  and  minimum  values  exhibit  “crankshaft”  chaos.  Without  precession,  the 
inclination  magnitude  is  much  smaller,  and  the  node  regresses  uniformly 


Goldreich’ s  approximation  does  not  adequately  model  the  actual  dynamics,  since  it  predicts 
a  much  larger  secular  change  than  the  semi-analytical  model. 

All  in  all,  the  dynamics  (i.e.,  particular  trajectories)  predicted  by  the  two  models  are 
quantitatively  different.  At  the  same  time,  when  it  comes  to  the  most  physically  important 
conclusion  from  the  Goldreich  approximation,  the  “Goldreich  lock”  of  the  inclination,  one 
may  still  say  that,  qualitatively,  the  semianalytical  model  confirms  the  locking  even  in  the  case 
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Fig.  10  Preparation  to  comparison  of  the  semianalytical  model  to  Goldreich’s  model.  In  Goldreich’s  model, 
the  Martian  equator  is  assumed  to  precess  uniformly,  while  the  semianalytical  model  takes  into  account 
(via  Colombo’s  equation)  variations  of  the  equinoctial  precession 


x  106 


Time  [yr] 


Fig.  11  Comparison  of  the  semi-analytical  model  to  Goldreich’s  model  for  a  satellite  of  Deimos’  mass  and 
with  the  initial  conditions  given  by  (61).  While  in  Goldreich’s  model  the  inclination  remains  tightly  locked, 
the  semianalytical  model  reveals  much  larger  variations  of  i 
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when  the  equinoctial  precession  is  variable,  and  the  solar  pull  on  the  satellite  is  included. 
The  locking  survives  even  for  highly  inclined  satellite  orbits.  It  is,  though,  not  as  stiff  as 
predicted  by  the  Goldreich  model:  we  can  see  from  Fig.  11  that  the  orbit  inclination  varies 
within  a  five-degree  span,  while  the  Goldreich  approximation  would  constrain  it  to  fractions 
of  a  degree. 

An  intriguing  fine  feature  of  the  inclination  evolution,  which  manifests  itself  for  orbits 
close  to  polar,  is  the  “crankshaft”.  This  kind  of  behaviour,  well  defined  in  Fig.  11,  is  not 
rendered  by  the  Goldreich  model,  because  that  model  was  initially  developed  for  small 
inclinations  and  uniform  precession.  One  might  suspect  that  the  “crankshaft”  is  merely  a 
numerical  artefact,  had  it  not  been  discovered  under  different  circumstances  (in  the  absence 
of  precession  but  in  the  presence  of  a  third  body)  by  Zhang  and  Hamilton  (2005).  This  kind 
of  pattern  may  be  generic  for  the  close  vicinity  of  i  =  90  deg. 


5  Conclusions 

In  the  article  thus  far,  we  continued  developing  a  tool  for  exploring  long-term  evolution  of  a 
satellite  orbit  about  a  precessing  oblate  primary.  In  particular,  we  were  interested  in  the  time- 
dependence  of  the  orbit  inclination  relative  to  the  moving  equator  of  date.  Our  model  includes 
three  factors:  J2  of  the  planet,  the  planet’s  nonuniform  equinoctial  precession  described  by 
the  Colombo  formalism,  and  the  gravitational  pull  exerted  by  the  Sun  on  the  satellite.  The 
problem  was  approached  using  different  methods.  One,  semianalytical,  was  based  on  numer¬ 
ical  integration  of  the  averaged  Lagrange-type  equations  for  the  secular  parts  of  the  Keplerian 
orbital  elements  introduced  in  a  noninertial  reference  frame  coprecessing  with  the  planetary 
equator  of  date.  The  right-hand  sides  of  these  equations  consisted  of  precession-generated 
contributions  and  contributions  due  to  the  direct  pull  of  the  Sun.  The  other  approach  was 
a  straightforward,  purely  numerical,  computation  of  the  satellite  dynamics  using  Cartesian 
coordinates  in  a  quasi-inertial  reference  frame. 

The  results  of  both  computations  have  demonstrated  a  good  agreement  over  10  Myr.  This 
means  that  the  semianalytical  model  adequately  describes  the  dynamics  over  this  time  span. 
Specifically,  the  terms  neglected  in  the  semianalytical  model  (the  short-period  terms  and  the 
terms  of  order  0(/l2))  play  no  significant  role  on  this  time  scale. 

Our  calculations  have  shown  the  advantages  and  the  limitations  of  a  simple  model  devel¬ 
oped  by  Goldreich  (1965)  for  uniform  equinoctial  precession  and  low  inclinations.  Though 
his  model  does  not  adequately  describe  the  dynamics  (that  turns  out  to  be  chaotic),  the  main 
physical  prediction  of  Goldreich’ s  model — the  “Goldreich  lock” — sustains  the  presence  of 
the  Sun  and  variations  of  equinoctial  precession,  provided  the  initial  inclination  is  suffi¬ 
ciently  low.  For  low  initial  inclinations,  the  inclination  exhibits  variations  of  order  fractions 
of  a  degree.  For  higher  inclinations,  however,  it  varies  already  within  a  span  of  about  ten 
degrees.  For  near-polar  orbits,  the  inclination  behaviour  demonstrates  the  “crankshaft”,  an 
chaotic  pattern  not  accounted  for  by  the  Goldreich  model.  The  “crankshaft”  emerges  because 
in  our  model  both  the  precession  variations  and  the  pull  of  the  Sun  are  included  into  the  model. 
However,  numerical  experiments  have  also  shown  that  the  “crankshaft”  gets  generated  by 
each  of  these  two  factors  separately. 
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Appendix  A:  Closed-form  expressions  for  /ti,  /62,  ii 3 

Using  the  compact  notation  cq  =  cos(-)  and  sq  =  sin(-),  and  conforming  to  the  procedure 
described  in  the  text,  we  obtain  the  following  expressions  for  /X2,  fiy. 


Ai  =  (((6  p2q2  ~  P4  ~  q4)(ci„f  +  ( -6p2q 2  +  p4  +  q4)cip)(chp)4 
+  ((-3  p4  -  6  p2q2  -  3  q2  +  3  p2  +  5  q4){ clp)3 
+  (2p4  -4  q4  +6  p2q2  -2p2  +lq2)cIp)  ■  (c (hp))2 
+  (3  q2  -  3  p2q2  -  4 q4)(cip)3  +  (3  q4  -  2  q2  +  2  p2q2) cIp 
+  (((4 qp3  -4q3p)(cIp)3  +  (4 q3p  -  4  qP3)cIp)(chp)5 
+  ((2  qp3  -6  qp+  14q3p)(cIp)3 

+  (~12q3p  +  4qp)clp)(chp)3  +  ((6 qp  -  10 q3 p  -6qp3)(cIp)3 
+  (4qp3  +  Sq3p  -4qp)cIp)chp)  ■  (shp)~l)(sIp)~l 

+  ((-9  pq2  +  9  pq4  -  3  p5  +  3  p3  +  6  p3q2)(cip )2  +  3  pq2  -  2  p3q2  +  p5 

-  3pq4  -  p3){ c3hp  +  ((-11  p3 q2  -  p5  -  10  pq4  -  p  +  2  p3  +  1 1  pq2)(c(Ip))2 

+  3pq4  +  3  p3q2  -  3pq2)chp  +  (((-9  p2q  +  6p2q3  +  9  p4q  -3  q5  +3q3)(clp)2 

-  3p4 q  +  3  p2 q  —  2  p2 q3 

~  q 3  +  q 5)(Ch„  +  ((7  P2q  -  P2q 3  -  8  q3  -  8  p4q  +  1  q5  +  q)(cip)2 

-  2q5  +  p2q 3  -\-2q 3 

+  3 p4q  -  3  p2q)(c(hp))2  +(5q3  -4  q5  -q  -5  p2q3  -  p4q  +  2p2q)(cIp)2 


-«  ^(((-2 qp  -2 q  p)( cIp)2  +  2 qp  +  2q  p)(chp)2 

+  (q  p  +  qp){cIp)2  -  qp  ~qp  +  (((-2  pp  +  2qq)(cIp)2  +  2  pp  -  2qq)(chp)3 
+  ((2  pp  -  2qq)(cIp)2  -2 pp  +  2qq)chp)(s(hp))~1)(s(Ip))~l 

+  ^(“4  p2  ~2q2q  +  q  -  qpp)cIpchp 

(A~2/>2A~Mg  -  Ag2)c/pfap)2  +  (~A  +  Ag2  +  2p2A  + 

Shp  ) 


(62) 


A2  =  y.(4q3p-4qp3)(cIp)2+4qp3 -4q3p)(chp)4 

+  ((qp2  -1  q3 p  +  2qp){cIp)2  -3qp3  +5q3p)(chp)2 
+  (2  q3p  —  qp  +  qP3)(cIp)2  -  q3 p  +  (((6  p2q2  -  p4  -  q4)(cIp )2 
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-  6 p2q2  +  p4  +  q4)(c5hp  +  ((3  q4  -  q1  -  9  p2q2  +  p2)( CIp)2 
~  P 4  -  2q4  +  9  p2q2)(chp)3  +  ((-p2  +  q2  -2  q4  +  p4  +  3  p2q2)(c(Ip))2 
+  q4  -  3  p2q2)c(hp))(shp)~l)(s(Ip))~l  +  ((-4p2q3  +  6 p2q  -2q3 
+  2q5  -6p4q)c(Ip)(chp)3  +  (2q3  +4  p4q  +  2p2q3  -4 p2q  -  2q5)cIpchp 
+  ((-2  p5  -  6  pq2  +  6  pq4  +  4  p3q2  +  2  p3)cIp (chp)4 
+  (2p5  +  8pq2-8  pq4  -  6  p3q2  -  2  p3)cIp  (chp)2  +  (2  pq4  -  2  pq2 


+  2p3q2)cIp)(shp)  b 


a 


2 


- a((2  pp  -  2  qq)(C[p)3  +  (-2  pp+  2 qq)cIp)(c2hp 
+  (2  pp  +  4qq)(C[p)3  +  (-2  pp  -  4qq)clp  +  (((-2qp  -  2q  p)(cIp)3 
+  (2qp  +  2q  p) cIp)(chp)3  +  ((2qp  +  2q  p)(c(Ip))3 
+  (-2 qp  -  2q  p)cIp)c(hp)){&hp)~l){&Ip)~x 

+  ((2  p  q2  +  2  pqq  +  4  p2  p  -  2  p)(cIp )2  +  p-  2  p2p  -  pqq  -  p  q2)  chp 
+  | \{~2qpp  -  2q  p2  -  4q2q  +  2q)(cIp)2  -  q  +q  p2  +  2q2q  +qpp)(c2hp 
+  (2qpp  ~  2q  +  2q  p2  +  4q2q)(cIp)2 


-  qp 


2 


2  q2q+q  -  qpp)  ■  (s (hp))  1 


(63) 


A 3  =  [(((4q3p-4qp3)(cIp)2  +4qp3  -4q3p)(chp)4  +  ((qp3  -  1  q3 p  +  2qp)(c(Ip))2 

-  3  qp3  +  5  q3p)(chp)2  +  (2  q3 p  -  qp  +  qp3)(cIp)2  -  q3 p 
+  (((6  p2q2  -  p4  -  q4)(cip)2  -  6 p2q2  +  p4 

+  q4)(c(hp))5  +  ((3  q4  -  q2  -9  p2q2  +  p2)(cIp)2  -  p4  -  2q4  +  9  p2q2)(chp)3 
+  ((-p2  +  q2  ~  2 q4  +  p4  +  3  p2q2)(c(Ip))2  +  q4  -  3  p2q2)c(hp))(shp)~1)(s(Ip))~l 
+((-4  p2q3  +  6  p2  q  —  2  q3  +2  q5 

-6p4q)c(Ip)(chp)3  +  (2q3  +  4 p4q  +  2p2q3  -4p2q  -2q5)cIpchp 
+  ((-2  p5  -  6  pq2  +  6  pq4  +  4  p3q2  +  2  p3) cIp (chp)4 
+  (2p5  +  8  pq2  -8  pq4  -6p3q2  -  2  p3)cIp(chp)2 

+(2  pq4-2  pq2+2  p3q2)cIp)(shp)~1)Jl-p2-q2  ^  a2-a  ^(((2  pp-2qq)(cIp)3 

+  (-2  pp  +  2qq)c(Ip))(c(hp))2  +  (2  pp  +  4qq)(cIp)3  +  (-2  pp  -4qq)cjp 
+  (((-2 qp  -  2qp)(cIp)3  +  (2 qp  +  2q  p)cIp)(chp)3  +  ((2qp  +  2q  p)(c(Ip))3 
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+  (-2  qp  -  2q  p)clp)c(hp))(shp)~l)(sIp)~l  +  (((2  pq1  +2  pqq  +  4p2p 

-  ip)(cip)2  +  p-  ip2p-pqq 

~  pq2)chp  +  (((-2  qpp - 2  q  p1  -4 q2q +2  q)(c(Ip))2 -q  +q  p2 +2  q2q 
+  qPP)(ch„)2+(2qpp-2q+2q  p2+4q2q)(cIp)2-q  p2 

-2q2q+q-qpp)(shp)~1)  ■  k/ 1  -  P2  ~  q2 


Appendix  B:  Niching  genetic  algorithms 

The  most  commonly  used  Genetic  Algorithm  (GA)  is  the  so-called  “Simple  GA”  (Goldberg 
1989).  To  perform  an  evolutionary  search,  the  Simple  GA  uses  the  operators  of  crossover, 
reproduction,  and  mutation.  A  crossover  is  used  to  create  new  solution  strings  (“children” 
or  “offspring”)  from  the  existing  strings  (“parents”).  Reproduction  copies  individual  strings 
according  to  the  objective  function  values.  Mutation  is  an  occasional  random  alteration  of 
the  value  of  a  string  position,  used  to  promote  diversity  of  solutions. 

Although  Simple  GA’s  are  capable  of  detecting  the  global  optimum,  they  suffer  from 
two  main  drawbacks.  First,  convergence  to  a  local  optimum  is  possible  due  to  the  effect  of 
premature  convergence,  where  all  individuals  in  a  population  become  nearly  identical  before 
the  optima  has  been  located.  Second,  convergence  to  a  single  optimum  does  not  reveal 
other  optima,  which  may  exhibit  attractive  features.  To  overcome  these  problems,  modifi¬ 
cations  of  Simple  GAs  were  considered.  These  modifications  are  called  niching  methods, 
and  are  aimed  at  promoting  a  diversity  of  solutions  for  multi-modal  optimisation  problems. 
In  other  words,  instead  of  converging  to  a  single  (possibly  local)  optimum,  niching  allows 
for  a  number  of  optimal  solutions  to  co-exist,  and  it  lets  the  designer  choose  the  appropri¬ 
ate  one.  The  niching  method  used  throughout  this  study  is  that  of  Deterministic  Crowding. 
According  to  this  method,  individuals  are  first  randomly  grouped  into  parent  pairs.  Each  pair 
generates  two  children  by  application  of  the  standard  genetic  operators.  Every  child  then 
competes  against  one  of  his  parents.  The  winner  of  the  competition  moves  on  to  the  next 
generation.  By  using  the  notation  P/  for  a  parent,  C;  for  a  child,  /(•)  for  a  fitness,  and  d(-) 
for  a  distance,  a  pseudo-code  for  the  two  possible  parent-child  tournaments  can  be  written 
as  follows  (Gurfil  and  Kasdin  2002a,  b): 

If  [d(Pu  Ci)  +  d(P2,  C2)  =  d{Pu  C2)  +  d(P2,  Ci)] 

If  /(Ci)  >  /(Pi)  replace  Pi  with  Ci 
If  f(C2)  >  f(P2 )  replace  P2  with  C2 
Else 

If  /(Ci)  >  /(P2)  replace  P2  with  Ci 
If  /(C2)  >  /(Pi)  replace  Pi  with  C2 

In  addition  to  applying  the  Deterministic  Crowding  niching  method,  we  used  a  two-point 
crossover  instead  of  a  single-point  one.  In  the  Simple  GA,  the  crossover  operator  breaks 
the  binary  string  of  parameters,  the  “chromosome,”  at  a  random  point  and  exchanges  the 
two  pieces  to  create  a  new  “chromosome.”  In  a  two-point  crossover,  the  “chromosome”  is 
represented  with  a  ring.  The  string  between  the  two-crossover  points  is  then  exchanged.  The 
two-point  crossover  or  other  multiple-point  crossover  schemes  have  preferable  properties 
when  optimisation  highly  nonlinear  functions  is  performed. 
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